Convolution


By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. 1D Convolution

Convolution is defined as the integral of the product of the two functions after one is reversed and shifted


$$y[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] = x[n]*h[n]$$

Output $y[n]$ came out by convolution of input $x[n]$ and system $h[n]$

  • Cross correlation and convolution
In [1]:
%%html
<center><iframe src="https://www.youtube.com/embed/Ma0YONjMZLI?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

2. 1D Convolution Examples

Convolution = filtering in time domain

  • We need to properly design the convolution filter to extract right features from inputs.

2.1. 1D Convolution

Installation guide

Type the follwoings in the cmd prompt

  • pip install ipython
  • pip install librosa
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import IPython.display as ipd
import librosa.display

from scipy import signal
In [3]:
N = 8
n = np.arange(N)
x = [0, 0, 0, 1, 1, 1, 1, 0]
h = [0, 0, 0, 0, 1/3, 2/3 ,1, 0]

# Convolve pulse with itself
y = np.convolve(x, h)

# plot
plt.figure(figsize=(10, 12))
plt.subplot(3,1,1)
plt.stem(x)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('x', fontsize = 15)
plt.subplot(3,1,2)
plt.stem(h)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('h', fontsize = 15)
plt.subplot(3,1,3)
plt.stem(y)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse * pulse', fontsize = 15)
plt.show()
In [4]:
# Convolve a rectangular pulse with itself

# pulse
N = 8
n = np.arange(N)
x = [0, 0, 0, 1, 1, 1, 0, 0]

# Convolve pulse with itself
y = np.convolve(x, x)

# plot
plt.figure(figsize=(10, 8))
plt.subplot(2,1,1)
plt.stem(x)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse', fontsize = 15)

plt.subplot(2,1,2)
plt.stem(y)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse * pulse', fontsize = 15)
plt.show()

Illustrative explanation


2.2. Denoising a Piecewise Smooth Signal

  • moving average (MA) filter

  • low-pass filter in time domain

moving average (MA) filter

  • a moving average is the unweighted mean of the previous $m$ data


$$\bar{x}[n] = \frac{x[n] + x [n-1] + \cdots + x[n-m+1] }{m}$$

In [5]:
# piecewise smooth signal
N = 50
n = np.arange(N)

v = np.hstack([np.ones([1, np.int(N/2)]), -np.ones([1, np.int(N/2)])])
x = np.sin(np.pi/N*n)*v
xn = x + 0.1*np.random.randn(1, N)

# construct moving average filter impulse response of length M
M = 5
h = np.zeros([1, M])
h[0:M] = 1/M

# convolve noisy signal with impulse response
y = signal.convolve(xn, h)

# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')
plt.subplot(2,2,2)
plt.stem(xn.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')
plt.subplot(2,2,3)
plt.stem(h.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')
plt.subplot(2,2,4)
plt.stem(y.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()

2.3. Edge Detection

  • Edge detection using first derivative


$$\frac{d}{d n} x[n] = x[n] - n[n-1]$$

In [6]:
# haar wavelet edge detector
M = 2
h = np.zeros([1, M])
h[0,0:1] = -1/M
h[0,1:2] = 1/M

# convolve noisy signal with impulse response
y = signal.convolve(xn, h)

# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')
plt.subplot(2,2,2)
plt.stem(xn.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')
plt.subplot(2,2,3)
plt.stem(h.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')
plt.subplot(2,2,4)
plt.stem(y.T, linefmt = 'b', markerfmt = 'w', basefmt = 'r')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()
In [7]:
%%html
<center><iframe src="https://www.youtube.com/embed/I3isHB9Iy7E?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
In [8]:
%%html
<center><iframe src="https://www.youtube.com/embed/ZW6pw89cuEA?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

2.4. Convolution on Audio

In [9]:
x, sr = librosa.load('./data_files/violin_origional.wav')
x = x/np.max(x)    # normalized

ipd.Audio('./data_files/violin_origional.wav', rate = sr) # play a wave file with sampling rate sr
Out[9]:
In [10]:
x_cor = x + 0.02*np.random.randn(*x.shape)

ipd.Audio(x_cor, rate = sr)
Out[10]:
In [11]:
# plot
plt.figure(figsize=(10,8))
plt.subplot(2,1,1)
librosa.display.waveplot(x, sr = sr)
plt.xlim([0, 5])
plt.title('Original')
plt.xlabel('')
plt.subplot(2,1,2)
librosa.display.waveplot(x_cor, sr = sr)
plt.xlim([0, 5])
plt.title('Corrupted')
plt.show()
In [12]:
print(x_cor.shape)
x_cor = x_cor[np.newaxis, :]
print(x_cor.shape)

(110250,)
(1, 110250)
  • Smoothing (low-pass filter)
In [13]:
M = 5
h = np.zeros([1, M])
h[0:M] = 1/M

y = signal.convolve(x_cor, h)
y = y/np.max(y)

ipd.Audio(y[0], rate = sr) # play a wave file with sampling rate sr
Out[13]:
In [14]:
# plot
plt.figure(figsize=(10,8))
plt.subplot(2,1,1)
librosa.display.waveplot(x_cor[0], sr = sr)
plt.xlim([0, 5])
plt.title('Original')
plt.xlabel('')
plt.subplot(2,1,2)
librosa.display.waveplot(y[0], sr=sr)
plt.xlim([0, 5])
plt.title('Convoluted')
plt.show()
  • Convolution with the impulse response
In [15]:
x, sr = librosa.load('./data_files/violin_origional.wav')
x = x/max(x)    # normalized

ipd.Audio('./data_files/violin_origional.wav', rate = sr)
Out[15]:
In [16]:
h, sr = librosa.load('./data_files/gunshot.wav')

ipd.Audio('./data_files/gunshot.wav', rate = sr)
Out[16]:
In [17]:
y = signal.convolve(x, h)
y = y/max(y)

ipd.Audio(y, rate = sr)
Out[17]:

3. Convolution on Images

3.1. Images

We wil show you what images look like in their raw form.



Installation guide

Install it via shell/command prompt:

  • pip install scikit-image
In [18]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import cv2
import skimage
from scipy import signal
In [19]:
# original image

img = plt.imread('.\image_files\lena_rgb.jpg')

plt.figure(figsize = (6,6))
plt.imshow(img)
plt.axis('off')
plt.show()
In [20]:
# image shape

img.shape
Out[20]:
(497, 497, 3)

Pictures are just another matrix. Although their shape may seem odd. Color pictures are a matrix with a shape of (pixel width $\times$ pixel height $\times$ 3)

You can think of this as a 3D matrix, or three matricies of (pixel width $\times$ pixel height) stacked on top of each other. The word 'Tensor' is often used to desribe these matricies generally.

In [21]:
red = np.zeros(img.shape, dtype = 'uint8')
red[:,:,0] = img[:,:,0]

plt.figure(figsize = (6,6))
plt.imshow(red)
plt.axis('off')
plt.show()
In [22]:
green = np.zeros(img.shape, dtype = 'uint8')
green[:,:,1] = img[:,:,1]

plt.figure(figsize = (6,6))
plt.imshow(green)
plt.axis('off')
plt.show()
In [23]:
blue = np.zeros(img.shape, dtype = 'uint8')
blue[:,:,2] = img[:,:,2]

plt.figure(figsize = (6,6))
plt.imshow(blue)
plt.axis('off')
plt.show()
In [24]:
# gray image

img = cv2.imread('.\image_files\lena_rgb.jpg', 0)
print(img.shape)

plt.figure(figsize = (6,6))
plt.imshow(img, 'gray')
plt.axis('off')
plt.show()

(497, 497)

3.2. 2D Convolution on Images

Convolution in 2D



Filter (or Mask, Kernel)

  • Modify or enhance an image by filtering
  • Filter images to emphasize certain features or remove other features

  • Filtering includes smoothing, sharpening and edge enhancement

  • Discrete convolution can be viewed as element-wise multiplication by a matrix




In [25]:
# noise image

img = cv2.imread('.\image_files\lena_sigma25.png', 0)
print(img.shape)

plt.figure(figsize = (6,6))
plt.imshow(img, 'gray')
plt.axis('off')
plt.show()

(512, 512)
In [26]:
# smoothing or noise reduction

M = np.ones([3,3])/9

img_conv = signal.convolve2d(img, M, 'same')
print(img_conv.shape)

plt.figure(figsize = (12,6))
plt.subplot(1,2,1)
plt.imshow(img, 'gray')
plt.title('Noisy Image')
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(img_conv, 'gray')
plt.title('Smoothed Image')
plt.axis('off')
plt.show()

(512, 512)
In [27]:
# original image

img = cv2.imread('.\image_files\lena.png', 0)
print(img.shape)

plt.figure(figsize = (6,6))
plt.imshow(img, 'gray')
plt.axis('off')
plt.show()

(512, 512)
In [28]:
# guess what kind of image will be produced after convolution

Mv  = np.array([[-1, 0, 1],
                [-2, 0, 2],
                [-1, 0, 1]])

img_conv = signal.convolve2d(img, Mv, 'same')

plt.figure(figsize = (12,6))
plt.imshow(img_conv, 'gray')
plt.title('Mv')
plt.axis('off')
plt.show()
In [29]:
# guess what kind of image will be produced after convolution

Mh = np.array([[-1, -2, -1],
               [0, 0, 0],
               [1, 2, 1]])

img_conv = signal.convolve2d(img, Mh, 'same')

plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
In [30]:
M = (Mv + Mh)/2

img_conv = signal.convolve2d(img, M, 'same')

plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
In [31]:
M = np.array([[0, -1, 0],
              [-1, 4, -1],
              [0, -1, 0]])

img_conv = signal.convolve2d(img, M, 'same')

plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()

Gaussian filter for burring images

In [32]:
# Gaussian filter
M = 1/273*np.array([[1,  4,  7,  4, 1],
                    [4, 16, 26, 16, 4],
                    [7, 26, 41, 26, 7],
                    [4, 16, 26, 16, 4],
                    [1,  4,  7,  4, 1]])

M = skimage.transform.resize(M, [15, 15], mode = 'reflect', anti_aliasing = True)
In [33]:
img_conv = signal.convolve2d(img, M, 'same')

plt.figure(figsize = (6,6))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
In [34]:
plt.figure(figsize = (14,6))
plt.subplot(1,3,1)
plt.imshow(img, 'gray')
plt.title('Input image')
plt.axis('off')
plt.subplot(1,3,2)
plt.imshow(M, 'gray')
plt.title('Image filter (15 x 15) ')
plt.axis('off')
plt.subplot(1,3,3)
plt.imshow(img_conv, 'gray')
plt.title('Convoluted image')
plt.axis('off')
plt.show()
In [35]:
img_rgb = cv2.imread('.\image_files\gala.jpg')
img_gray = cv2.imread('.\image_files\gala.jpg', 0)
print(img.shape)

plt.figure(figsize = (18,10))
plt.imshow(img_rgb)
plt.axis('off')
plt.show()

(512, 512)
In [36]:
plt.figure(figsize = (18,10))
plt.imshow(img_gray, 'gray')
plt.axis('off')
plt.show()
In [37]:
# Gaussian filter
M = 1/273*np.array([[1,  4,  7,  4, 1],
                    [4, 16, 26, 16, 4],
                    [7, 26, 41, 26, 7],
                    [4, 16, 26, 16, 4],
                    [1,  4,  7,  4, 1]])

M = skimage.transform.resize(M, [23, 23], mode = 'reflect', anti_aliasing = True)

img_conv = signal.convolve2d(img_gray, M, 'same')

plt.figure(figsize = (18,10))
plt.imshow(img_conv, 'gray')
plt.axis('off')
plt.show()
In [38]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')